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This paper describes the implementations of the linear-eddy model (LEM) and an Eulerian 
FDF/PDF model in the National Combustion Code (NCC) for the simulation of turbulent 
combustion. The impacts of these two models, along with the so called ‘laminar chemistry’ 
model, are then illustrated via the preliminary results from two combustion systems: a nine- 
element gas fueled combustor and a single-element liquid fueled combustor. 


Introduction 

Recent more stringent emission standards have enticed the development of more fuel-efficient and low-emission 
combustion system for aircraft gas turbine applications. The analysis and design of such advanced combustors is 
facilitated by employing higher fidelity, hence, more sophisticated turbulent combustion models. With advances in 
computational power and the availability of distributed computers, the use of complex turbulence-chemistry 
interaction models is now becoming more and more affordable for realistic geometries. 

It is well known that the major difficulty, when modeling the turbulence-chemistry interaction, lies in the high 
non-linearity of the reaction rate expressed in terms of the temperature and species mass fractions. The so called 
‘laminar chemistry’ model uses directly the averaged temperature and averaged species mass fractions to calculate 
the averaged reaction rate (which is, at the best, valid only for low turbulence flows). The transport probability 
density function 1 (PDF) model and the linear eddy model 2 (LEM), both use local instantaneous values of the 
temperature and mass fractions, have been shown to often provide more accurate results of turbulent combustion. A 
filtered density function (FDF) is the spatial average of a PDF over a grid-box-sized volume 3 . The FDF has all the 
properties of a PDF 4 ; i.e., FDF is non-negative everywhere and is normalized. The transport equation for the 
FDF/PDF and others like it has been widely studied by combustion engineers since the pioneering work of 
Lundgren 5, 6 . The linear eddy model (LEM) represents an alternate strategy in turbulent combustion. It is a one- 
dimensional mixing model that combines molecular processes in a deterministic way (through the solution of 
unsteady reaction-diffusion equations), and turbulent transport through stochastic stirring events. 

The main purposes of this report are to describe the practical implementations of the LEM and the FDF/PDF 
models, to illustrate the differences among the LEM, the transport FDF/PDF and the ‘laminar chemistry’ via the 
simulation results of a hydrogen gas fueled combustor, and to demonstrate their applications to a Jet-A liquid fueled 
combustor. These preliminary results will then lead to further, more quantitative investigations to be reported in the 
near future. 


Governing Equations for Gas-phase Flow 

The governing equations for the (averaged) flow in the multiple rotating frames of reference (MRF), which are 
convenient for including the rotational components of the jet engine in the computational domain 7 , are summarized 
in a hybrid formulation below. 

Continuity equation: 


1 Aerospace Engineer, Combustion Branch, Member AIAA. 

2 Aerospace Engineer, Combustion Branch, Associate Fellow AIAA. 

1 

American Institute of Aeronautics and Astronautics 



j,\\\s’ dv+ \\f^-°,> dA =\\\ r s ^ dv <» 

Species transport equation: 
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V is an arbitrary control volume with control surface A , and all the flow variables as well as the source terms 
including the reaction rates are in some averaged sense, p is the fluid density, p m is the partial density of species 

m , w m is the production rate of species m per unit volume, u is the flow velocity in a stationary Cartesian 

coordinate system, E) g is the grid rotating velocity for each designated axis, U g = cb g x r , E is the total energy 

density, k is the turbulent kinetic energy, £ is the turbulent dissipation rate, S k is the source term for k , S £ is 

the source term for £ , p is the static pressure, T is the effective stress tensor, q is the effective heat flux vector, 

D eff is the effective mass diffusivity coefficient which is equal to EltL , p and jil t are the laminar viscosity and 

turbulent viscosity, respectively. S mlc , S mls , S mlm and S mle are source terms due to the contribution of spray, 

respectively. For gas, the Schmidt number (Sc ) is about 1. Different angular velocities (and even different rotating 
axes) are assigned to different mesh blocks or groups within the model. Balance equations for each group are 
expressed in the relative reference frame but in terms of the absolute velocity (i.e. the velocity with respect to a 
stationary coordinate system). 


Turbulence-Chemistry Interaction 


A rudimentary model for the turbulence-chemistry interaction is the so called ‘laminar chemistry’ model, which 
directly uses the averaged temperature and averaged species mass fraction to calculate the averaged reaction rates. 
Due to the high non-linearity of the reaction rates with respect to the temperature and species mass fractions, the 
laminar chemistry model is, at the best, valid only for low turbulence flows. This motivates the use of much more 
sophisticated models such as the FDF/PDF model and the LEM model. 

Linear-Eddy Model 

The linear-eddy model was proposed by Kerstein 8 for non-reacting flow and extended to reactive flow 9 . It has 
been investigated in detail in e.g. Reference [2], therefore, only a briefly summary is given here. The LEM approach 
accounts for both deterministic and stochastic contributions. The effects of the unresolved part of the velocity on the 
momentum and scalar transport are explicitly represented by turbulent stirring events. The LEM model is 
implemented in terms of a fractional splitting technique; it is divided into two processes: subgrid and supergrid. The 
subgrid process, which occurs within each computational cell, composes of four operators: (a) molecular diffusion, 
(b) finite-rate kinetics, (c) volumetric heat release, and (d) stochastic stirring that actually is a turbulent convection 
by a unresolved velocity component. The supergrid process represents the convection of the scalar field by the 
resolved velocity field. The subgrid process is achieved by first subdividing each cell into a fixed number of lumps, 
then the lumps are aligned from inflow to outflow within a cell forming a one dimensional path. The subgrid LEM 
governing equations are solved on each LEM lump. 

Subgrid processes: 

Species transport equation for LEM model: 
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F t stirring and F Y stirring are the unresolved turbulent convective terms and implemented explicitly using the 

stochastic re-arrangement events called triplet maps. The triplet map mimics the characteristics of the stretching and 
folding processes inherent to turbulent flows. The effect of a triplet map is a three-fold compression of the scalar 
field in a typical eddy size, / . The map is parameterized by the length-scale of a notional “eddy” acting on the 
lumps, t , the starting position in the lumps and the stirring frequency, X . For convenience, the position is randomly 
taken to indicate the starting location of the eddy from a uniform distribution. The eddy size is randomly selected 
from a prescribed distribution function given by 8 


/(/) = (5/3)r 8/3 /[/7“ 5/3 - C 5/3 ] 


( 9 ) 


Here the Kolmogorov scale 77 is determined from inertial range scaling law, 7j ~ l t Re 3/4 , l t is turbulence length 
scale. The event rate (frequency per unit length) is determined by 8 



The model constant C A is 15. The stirring time interval between events is give as a t stir = l/(/lA) where A is the 
total length of the lumps which is volume 173 of a cell. 

Supergrid processes: 

Lagrangian convection (also known as splicing) 


= - V -(u-jjg) (ll) 

dt g 

where (j) is Y { or T . This equation describes the mesh-resolved convection of the scalar field, and is solved via a 
Lagrangian transfer of mass across the finite-volume cell surfaces 10, n . 

Referring to Figure 1, for example, the outward mass through the right side of a cell computed from resolved 
velocity and density is equivalent to 1.5 LEM lumps and colored in red. The outward mass through the bottom side 
of the mesh is equivalent to 2.5 LEM lumps and colored in magenta. Similarly the inward mass through the top and 
left sides of the mesh are equivalent to 6 LEM lumps that are added to existing lumps. The net number of lumps is 
14 after the splicing step. It leads to a non-uniform number of lumps in each mesh. To avoid programming 
complexities in a parallel environment, the LEM domain is re-gridded to have lumps of equal volume and equal 
number after each splicing. 
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Arrows -- Convection direction 
Initial number of lumps - 12 
Light Green - inward mass (1 .5) 

Dark Green - inward mass (4.5) 

Red - outward mass (1 .5) 

Magenta - outward mass (2.5) 

Current number of lumps -- 14 

Final number of lumps will became 12 after re-gridding 
via equalizing the volume. 


Figure 1 Schematic illustrating the splicing algorithm used for scalar convection. (Assuming number of 
LEM lumps is 12. Each lump contains volume and species density) 


Some salient features of the LEM module are summarized below: 

1 . The LEM code receives the mean velocity, pressure and turbulence fields from the flow solver. And, also, 
it receives the source terms arising from the liquid-phase contribution from the spray solver if needed. 

2. The LEM solver provides the species and temperature solution to the flow solver and, also, to the spray 
solver if needed. 

3. The LEM model treats the nonlinear reaction rates without any approximation. 

4. The diffusion term in the subgrid process puts constraints on time step and usually requires multiple sub- 
iterations. 


FDF/PDF model 

The Eulerian joint scalar filtered mass-density function (FDF) approach for unsteady flow or the Eulerian joint 
scalar probability density function (PDF) for steady flow has been discussed in detail in the literature 1, 12, 13 and , 
therefore is only briefly summarized here. Like LEM mode, the transport equation of the probability density 
function is solved by making use of an approximate factorization scheme. It is split into four operators associated 
with chemical reactions, molecular mixing, spray and spatial transport. It is solved by making use of a Monte Carlo 
technique. In the Monte Carlo simulation the density function is represented by a large ensemble of stochastic 
particles. Each particle carries enthalpy, temperature and species mass fraction 

Subgrid processes: 

Species transport equation for FDF/PDF model 


^- = ^—C t aKr-Y,) + S„„ 
at p f 

Energy equation 

^ = -C t a,(h,- h) + S mk 


( 12 ) 


(13) 


Where CO = S Ik , C ^ is an empirical constant with a value of about one. The first term on the right hand side of the 
energy equation represents molecular mixing which is accounted for by making use of the relaxation to the ensemble 
mean submodel 14 . The last term, S mle , represents the contribution from the spray source terms. 

Supergid processes: 
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Eulerian convection: 


dl 

8t 


p Sc, 


( 14 ) 


where (/) is Y t or h £ , jU e is effective viscosity, Sc ^ is Schmidt number of (j) , p is the mixture density in the cell. 


Referring to Figure 2, unlike the LEM splicing algorithm, the Eulerian convection and diffusion process is achieved 
through the content replacement of particles, as each Monte Carlo particle carries the information of the mass 
fraction of species, the enthalpy, the temperature, but not the volume. For example, the net inward mass through the 
right side of a cell, computed from upwind-resolved convective terms and central-differenced effective diffusion 
term; then normalized by the mass in the cell, is equal to 5% and colored in red. The net inward mass through the 
bottom side of the mesh is equal to 1 0% and colored in magenta. Similarly the inward mass through the top and left 
sides of the mesh are equivalent to 40% and colored in dark and light green. If the total number of the particles is 20 
in this cell, then one particle (5% of 20) will be randomly selected from its right adjacent cell, and its contents are 
copied to one randomly selected particle in the current cell. The same procedure is applied to other particles in 
colors. For those remaining particles (45% here) the contents are randomly shuffled within the cell. This is 
equivalent to an upwind differencing scheme. 



Arrows -- Convection direction 
Light Green - inward mass (15%=3) 

Dark Green - inward mass (25%=5) 

Red - net inward mass (5%=1 ) 

Magenta - net inward mass (10%=2) 

Number of Monte Carlo particles - 20 
Contents of randomly selected particles are 
replaced by those of randomly selected particles 
from adjacent cells. 


Figure 2. Schematic illustrating the process of the scalar convection and diffusion for the Eulerian scalar 
Monte Carlo FDF/PDF model. 


Some salient features of the FDF/PDF module are summarized below: 

1 . The FDF/PDF code receives the mean velocity, pressure and turbulence fields from the flow solver. And, 
also, it receives the source terms arising from the liquid-phase contribution from the spray solver if needed. 

2. The FDF/PDF solver provides the species and temperature solution to the flow solver and, also, to the spray 
solver if needed. 

3. The FDF/PDF model treats the nonlinear reaction rates without any approximation. 


Finally, all the operators of LEM or FDF/PDF are integrated by an explicit Euler method. 


Preliminary Results 


The effects of employing different turbulence-chemistry models for turbulent combustion are illustrated by the 
calculated results from the simulations of two combustors: a Sandia 9-element hydrogen combustor and a NASA 
single-element lean direct injection combustor. 
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Sandia 9-element hydrogen gas fueled combustor 


The simulation of the hydrogen burner involved two separate computations. First a simulation of just one nozzle 
was performed, which included an upstream settling chamber, and a downstream expansion. The size of the settling 
chamber and expansion chamber was made to provide an identical expansion of area as that of the nine-nozzle case. 
The two hydrogen jets were provided with plenums, to damp out any possible fluctuations of the inlet conditions. 

After the simulation of the single nozzle case had run to convergence, data was taken at a location .005 m. 
upstream of the expansion region. This data is just downstream of the hydrogen jets, and is used as the initial 
conditions for the nine-nozzle simulation. 


The quantities recorded at this location are the distributions of the velocity, turbulence and scalar variables over 
a 2D area. These fields are then interpolated onto the inlet regions of the nine-element simulation. Each of the nine 
nozzles had the initial data so orientated that the hydrogen jets corresponded to the experimental location. The grid 
used contained 282,624 hexahedrons and is shown in Fig. 3. The chemical kinetics is represented by a 9-species, 20- 
step mechanism for hydrogen-air combustion. 


Z-O OC4 plane 


Figure 3. Schematic of the mesh used in the simulation of a Sandia 9-element hydrogen burner. It consists of 
282624 hexahedrons. Length is about 0.05 meter and diameter is about 0.038 meter. Each injector is 0.076 
meter in diameter. 

Computed axial velocity on the vertical center plane (x=0.0 plane) by using different turbulent combustion models 
are presented in Figure 4. The top two plots are results of URANS simulations while the bottom three plots are 
results of RANS simulations. It can be seen that the pattern of axial velocity from URANS and RANS is quite 
different. Figure 5 shows the axial velocity on z=0.004 m plane viewed from the exit. It indicates that results 
between FDF/PDF and laminar chemistry are quite different too. Figure 6 presents the pressure contours on x=- 
0.0014m. 

The temperature distribution on z=0.004 m is depicted in Figure 7. The results from the laminar chemistry model are 
distinctly different from those from the FDF/PDF model. Figure 8 shows that RANS simulations yield lower 
temperature flows from the inlet of the central injector all the way to the exit of the combustion chamber, but this is 
not the case for URANS simulations. The mass fraction of OH, a minor species, is presented in Figure 9 (z=0.004 m 
plane) and Figure 10 (x=0.0 m plane). Basically, it reflects the structure of the temperature. 
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Figure 4. Axial velocity contours at x=0.0 m plane. Plots in clockwise direction: FDF/URANS with steady 
pressure exit; FDF/URANS with unsteady pressure exit; PDF/RANS with sensible enthalpy correction; 
PDF/RANS; and laminar chemistry/RANS. 


axial velocity, M/S 
k 460 




Figure 5. Axial velocity contours at z=0.004 m plane. Plots in clockwise direction: FDF/URANS with steady 
pressure exit; FDF/URANS with unsteady pressure exit; PDF/RANS with sensible enthalpy correction; 
PDF/RANS; and laminar chemistry/RANS. 
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Figure 6. Pressure contours at x=0.0014 m plane. Plots in clockwise direction: FDF/URANS with steady 
pressure exit; FDF/URANS with unsteady pressure exit; PDF/RANS with sensible enthalpy correction; 
PDF/RANS; and laminar chemistry/RANS. 



Figure 7. Temperature contours at z=0.004 m plane. Plots in clockwise direction: FDF/URANS with steady 
pressure exit; FDF/URANS with unsteady pressure exit; PDF/RANS with sensible enthalpy correction; 
PDF/RANS; and laminar chemistry/RANS. 
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Figure 8. Temperature contours at x=0.0 m plane. Plots in clockwise direction: FDF/URANS with steady 
pressure exit; FDF/URANS with unsteady pressure exit; PDF/RANS with sensible enthalpy correction; 
PDF/RANS; and laminar chemistry/RANS. 
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Figure 9. Species OH mass fraction contours at z=0.004 m plane. Plots in clockwise direction: FDF/URANS 
with steady pressure exit; FDF/URANS with unsteady pressure exit; PDF/RANS with sensible enthalpy 
correction; PDF/RANS; and laminar chemistry/RANS. 
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Figure 10. Species OH mass fraction contours at x=0.0 m plane. Plots in clockwise direction: FDF/URANS 
with steady pressure exit; FDF/URANS with unsteady pressure exit; PDF/RANS with sensible enthalpy 
correction; PDF/RANS; and laminar chemistry/RANS. 


So far, we have presented the results obtained in the context of RANS and URANS simulations. In the following, 
the results obtained in the context of the partially resolved numerical simulation for a very large eddy simulation 
(PRNS/VLES) 15, 16, 17 will be presented. The PRNS/VLES is based on the concept of temporal filtering, in which 
the contents of both the resolved and unresolved turbulence are regulated by a “resolution control parameter” (RCP), 
which is related to the width of the temporal filter. A value of RCP =0 .4 is used in the following cases. 

While the PRNS is used to solve the pressure (i.e. the continuity equation), the momentum and the turbulence 
equations, the temperature (or enthalpy) and the species mass fraction are solved either via the FDF or the linear 
eddy model (LEM) approach. For sub-grid scalar field evolution, 24 LEM cells are used within each PRNS cell, 
allowing a resolution in the sub-grid approximately of the order of 77 , the Kolmogorov scale. For the scalar FDF 
simulation, 100 particles were used in each PRNS cell to represent the joint scalar FDF. 

In Figure 11, several isosurface plots of the axial velocity computed by LEM/PRNS are shown. It indicates that the 
center injector has the largest penetration. The distribution of the OH mass fraction is depicted in Figure 12. 
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Figure 11. Isosurfaces of the axial velocity (partial view). LEM/PRNS with unsteady pressure exit. 
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Figure 12. Isosurfaces of the OH mass fraction (partial view). LEM/PRNS with unsteady pressure exit. 
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Selected results from the FDF/PRNS and the LEM/PRNS calculations are shown in Figure 13 to 15. Figure 13 
shows the isosurface of axial velocity at 200 m/s by FDF (top plot) and by LEM (bottom plot) approaches. It 
suggests that the flow field from LEM is more organized than that from FDF. Figure 14 depicts the isosurface of 
temperature at 2000K. It suggests that the penetration of each injector by FDF approach is much deeper than that of 
LEM. It is speculated that LEM might provide stronger mixing. In Figure 15, the mass fraction isosurface of the 
species OH is shown, and it resembles the temperature structure. 
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Figure 13. Axial velocity isosurface at the value of 200 m/s. Top: FDF/PRNS with unsteady pressure exit; 
Bottom: LEM/PRNS with unsteady pressure exit. 
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Figure 14. Temperature isosurface at a value of 2000 K. Top: FDF/PRNS with unsteady pressure exit; 
Bottom: LEM/PRNS with unsteady pressure exit. 
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Figure 15. OH mass fraction (.004) isosurface. Top: FDF/PRNS with unsteady pressure exit; 
LEM/PRNS with unsteady pressure exit. 


Bottom: 
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NASA single-element Jet-A liquid fueled LDI combustor 


Simulations of a single-element lean direct injection (LDI) combustor with liquid spray have been performed. 
Detailed description of the spray equations and solver can be found in a previous report 18 . The description of the 
geometry and computational grid has been reported in the literature 19 . Only the boundary conditions for the fuel 
injector are summarized here. The liquid fuel, Ci 2 H 2 3 , is injected into the domain from a pressure swirl injector. The 
diameter of the orifice is .0006 m. The fuel atomizer is set at a pressure of 1 lOkPa and delivers a flow rate of 0.025 
kg/min. The spray has a 90 degree spray angle which is the same as the converging-diverging venture attached to the 
swirler. 

In the laminar chemistry/PRNS case, an inflow droplet size distribution 20 is prescribed, 


dn 

n 


l3.5 


4.21x10 [ ]-e 

d 


— 16 * 98 (— )°' 4 dd 


32 


*32 


where n is the total number of the droplets and dn is the number of droplets in the size range between d and 
d + dd. This correlation also requires the specification of a Sauter mean diameter, d 32 , and the number of droplet 
classes. These specified inflow droplets will undergo evaporation without secondary breakup. 

In the FDF/PRNS case, instead of specifying the inflow droplets, a sheet breakup model 21 is used for the primary 
atomization; furthermore, an ETAB (enhanced TAB) model 22 is used for the secondary droplet breakup. The value 
of RCP is set to 0.4 whose selection was discussed in Reference 23. It is also noted here that the solution of the 
laminar chemistry/PRNS simulation is used as the starting condition for the FDF/PRNS simulation. 

Snap shots of instantaneous axial velocity contour at the center plane from PRNS simulations of LDI are shown in 
Figure 16. The variation of the axial velocity in the case of using FDF is more noticeable than the variation in the 
case of using laminar chemistry. The temperature results are shown in Figure 17. The irregularity is also more 
pronounced in the case of using FDF. In Figurel8, the sliced contours of the instantaneous temperature at 0.04 meter 
downstream of the exit of the fuel injector are overlaid by instantaneous temperature isosurfaces at the same 
location. It can be seen that the results from using the FDF model indicate a more irregular swirling and an earlier 
merging of the lower temperature portion of the flow. 

In the LEM/PRNS case, due to its robustness the prescribed particle size distribution is used as the inflow boundary 
condition for the fuel injector. For sub-grid scalar field evolution, 18 LEM cells are used within each PRNS cell and 
up to 40 sub-iteration steps are set for the stirring event and molecular diffusion. The isosurfaces of the mean axial 
velocity in Figure 19 shows the zone of the recirculation has rather strong reverse flow, about 20 m/s. In Figure 20 
the mean turbulence kinetic energy is mild at the axis near the injector and reduces to smaller values as the flow is 
marching down stream. The mean temperature at the center plane and several axial stations are shown in Figure 21 
and 22 respectively. The influence of the swirler is clearly seen even at the exit plane. In the recirculation zone, the 
mean temperature at the axis is higher than that of surrounding regions. The equivalence ratio computed from the 
inlet boundary conditions is about 0.72. The adiabatic flame temperature is around 2100 K. The mean temperature at 
the exit plane is between 1800 K and 1900 K. 
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Axial velocity, m/s 



Figure 16. Snap shots of instantaneous axial velocity contour at the center plane from PRNS simulations of 
LDI: (top) laminar chemistry, (bottom) FDF. 
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Figure 17. Snap shots of instantaneous temperature contour at the center plane from PRNS simulations of 
LDI: (top) laminar chemistry, (bottom) FDF. 
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Figure 18 . Instantaneous temperature isosurface and sliced contours at -0.04 m downstream of the fuel 
injector (viewed from injector toward exit): left plot from laminar chemistry/PRNS, right plot from 
FDF/PRNS. 
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Figure 20. Mean turbulence kinetic energy isosurface at the center plane from PRNS simulations of LEM. 



Figure 22. Mean temperature contour plots at six axial stations from PRNS simulations of LEM. 
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Concluding Remarks 


Higher fidelity models for turbulence-chemistry interaction are implemented into the NCC. The differences 
among the LEM, the transport FDF/PDF and the ‘laminar chemistry’ are examined using the simulation results of a 
hydrogen gas fueled combustor. Their applications to a Jet-A liquid fueled combustor are also demonstrated. These 
preliminary results provide the basis for further, quantitative investigations of these higher fidelity models when 
applied to two-phase turbulent combustion systems. In Part II of the assessment task, the results of the EDI 
simulations from different turbulent combustion models will be compared with available experimental data. 
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